Taxonomic profiling

We used Kraken2/Bracken to perform taxonomic assignment of the reads. The procedure included the remotion of reads mapping the rice genome (O. zativa, Build_4.0). The workflow is proposed by Lu et al. (2022) as follows (wrapped by Nextflow):

The analysis of the raw counts was developed by wrapping the output as a Phyloseq object to estimate effect of treatments and time over relative abundance, α-Diversity, β-Diversity and differentially abundant taxa.

Metagenome-Assembled Genome (MAG) recovery

nf-core/mag was implemented to build and annotate the MAGs from the raw sequences. A co-assembly/co-binning strategy was chosen as an attempt to improve the quality and resolution of the MAGs. MEGAHIT was used for the assembly, while MaxBin2 and MetaBAT2 were the selected tools for binning. A de-replication and refinement step was carried by enabling DASTool within the workflow. For the downstream analysis, the curated MAGs were analyzed with MAGFlow/BIgMAG, that includes quality-measuring tools such as CheckM2, BUSCO, GUNC and QUAST, as well as the taxonomic annotator GTDB-Tk2 relying on the database version R220. The presence of carbohydrate-related metabolism enzymes was studied through dbCAN3 along with the CAZy database, whilst Cluster of Orthologous Genes (COGs) were detected by COGclassifier. The selection of the MAGs was carried out using the criteria: completeness above 50% and contamination below 10%. In addition, MAGS belonging to either rice straw or soil with average differences <5% ANI were clustered using FastANI. For data visualization, Circos, iTOL and ggplot2 packages were used.

Relative abundance, α-Diversity and β-Diversity

Figure 1. Taxonomic profiling of the sequences using Kraken2/Bracken: a Relative domain abundance per matrix source, merging all the replicates per sample. b α-Diversity, colored by Treatment, measured through Shannon and Inversed-Simpson indices. c Sample ordination through Principal Coordinate Analysis (PCoA) using a Bray-Curtis dissimilarity matrix as input. The ellipses represent 95% confidence regions around the centroids of each group calculated by the assumption of multivariate t-distribution. d Two-dimensional Uniform Manifold Approximation and Projection (UMAP) embedding of β-Diversity estimation of the samples based on the Aitchison distance. e Compositional variability of the samples, measured as distance to the centroid in a multivariate space computed with the Aitchison distance among samples. ANOVA analysis results show no significant differences in terms of group dispersion. f Proportion of the variance explained by the factors Time and Treatment, segmented by matrix source, in Permutational Analysis of Variance (PERMANOVA) implemented according to different distance/dissimilarity metrics. The significance of the each factor is determined by the symbol color. All the analysis from b to f were carried out only with Bacteria data and with normalized raw counts as Counts per million (CPM). The Treatment colors are transversal to plots b, d, e and f.

Figure 2. Relative phylum (a, c) and genus (b, d) abundance of agglomerated data for Bacteria (top) and Eukaryota (bottom). The displayed abundance is part only of the portion of the reads belonging to each domain ( see Figure 1a and Figure S1). Replicates were merged, and the panels are divided by matrix source. On phylum plots, only phyla with relative abundance above 0.5% are displayed, whilst on genus plots genera above 3% and 2% are depicted for Bacteria and Eukaryota, respectively.

Differentially abundant taxa

Figure 3. Selected differentially abundant genera in rice straw (a) and soil (b) samples based on Analysis of Compositions of Microbiomes with Bias Correction 2 (ANCOM-BC2). Right: Log2FC of Treatments to Control. Left: Boxplots of selected genera abundance (counts) in log10 scale. All displayed genera exhibit p-value < 0.05 adjusted with Holm–Bonferroni method for multiple comparisons. Cont0 and Cont1 in rice straw comparisons indicate Control samples at Time 0 and 1, respectively. As a result, rice straw results depict differences throughout time for all treatments including Control. In the case of soil samples, Time 0 and Time 1 were merged for each Treatment since Permutational Analysis of Variance (PERMANOVA) results does not show significant differences for this factor. ANCOM-BC2 was implemented only with Bacteria data, agglomerated at genus level and taxa abundant in at least 0.5% and 0.05% (any sample) for rice straw and soil samples, respectively.

Summary and presence/absence of genera across samples according to GTDB-Tk2 annotation

Figure 4. Metrics summaries of the MAGs built from rice straw (a) and soil (c) samples; the columns represent, in order, the proportion of: mid-quality and high-quality MAGs according to CheckM2 estimations; MAGs passing the GUNC test; and annotated MAGs at different taxonomical levels based on GTDB-Tk2. Presence/absence of MAGs annotated with GTDB-Tk2 in rice straw (b) soil (d) samples at genus level. Genome size, in Mbp, is indicated as the point size, while the color scale represents the relative abundance of each MAG per sample, and the point shape shows whether the displayed MAG is classified as high-quality or mid-quality (CheckM2). When MAGs with the same genus annotation were found in the same sample, the most complete one is depicted and their abundances are added up.

Cross-referencing of quality metrics, taxonomical annotation, COG estimation and CAZy information for clustered MAGs.

Figure 5. Recovered MAGs from rice straw (a) and soil (b) samples after clustering based on Average Nucleotide Identity (ANI). The visualization includes phylogenetic relationship among them (tree), specific highlighted phylum clades (colored sectors), deepest annotated taxon per bin (tips), genus classification (outer text), Treatment (color strip), Time information (black triangle binary), proportion of the main categories of Cluster of Orthologous Genes (COG) within each MAG (blue-orange-red-green pie chart), Carbon Active EnZyme (CAZy) class normalized abundance (red scaled heatmap), classification as either high-quality or mid-quality MAG based on checkM2 measurements (blue binary), GUNC test passing or not (green binary), and genome size in Mbp (outer black bars).The displayed MAGs are representatives from the generated clusters using ANI, and they were selected based on the following score: completeness – 0.5 x contamination. Conventions in CAZy heatmap: GH glycoside hydrolases, GT glyco-syltransferases, PL polysaccharide lyases, CE carbohydrate esterases, CBMs carbohydrate-binding modules, AA auxiliary activities.

MAG clustering using CAZy information and β-glucosidase activity estimation

Figure 6. a Two-dimensional Uniform Manifold Approximation and Projection (UMAP) embedding of normalized CAZy enzymes counts per MAG (before clustering with FastANI) colored by phylum (GTDB-Tk2). All CAZy classes (GH, GT, PL, CE, AA and CBMs), along with all their families, were considered for clustering using default values (15 neighbors and Euclidean distance). b Left: bubble plot depicting normalized CAZy counts for β-glucosidase families (GH1, GH3, GH5 and GH9) per sample; GH glycoside hydrolases. Right: data distribution for experimental β-glucosidase activity quantified per sample. The point size on the bubble plot represents the maximum value of the counts of each CAZy family per sample, and Treatment colors are shared by both the bubble plot and the boxplot.

Abundance of CAZy-derived lignocellulose-degrading enzymes per MAG

Figure 7. Bubble plot depicting the number of normalized counts for CAZy families associated with lignocellulose-degrading enzymes found in the MAGs recovered from rice straw (top) and soil (bottom) samples. Visualization includes (from left to right): MAG clustering using the CAZy counts of the selected families, Treatment origin of the MAG, phylum information (GTDB-Tk2), Time (0: initial, 1: after 30 days of Treatment). The classification as mid-quality or high-quality MAG (CheckM2) is represented by the shape in the bubble plot. The displayed MAGs are representatives from the generated clusters using ANI, and they were selected based on the following score: completeness – 0.5 x contamination. Additionally, for rice straw, MAGs were filtered considering only the ones with at least 3000 normalized counts in 5 or more CAZy families; complete information about normalized counts of these MAGs is found in figure S4.

Percentage of annotated reads with Kraken2/Bracken

Figure S1. Percentage of annotated reads using Kraken2/Bracken segmented by matrix and time point. For rice straw at Time 0 samples were taken from a bulk of material; this same bulk was subsampled with specific amounts of rice straw to start the experiments, and hence accounting with only 1 sample at Time 0.

Quality control with mock community

Figure S2. Krona plot of displaying the taxonomical classification and abundance (Kraken2/Bracken) of the reads belonging to the **mock community* sequenced as part of quality control for the DNA extraction process.

GTDB taxonomy of rice straw MAGs

GTDB taxonomy of soil MAGs

Cluster of Orthologous Genes (COG) annotation

Figure S3. Circos plots showing the proportion of the 4 main categories (Metabolism; Information Storage and Processing; Cellular Processes and Signaling; and Protein Coding Sequences with Unknown function) of Clusters of Orthologs Genes (COGs) found within the MAGs recovered from rice straw (a) and soil (b) the samples. The inner ring accounts for the total number of COGs related to every sample or COG category; the outer ring depicts the relative abundance of COGs from each sample or COG category; the width of the ribbons linking any sample and a COG category depicts their relative abundance to each other. The number of counts per category was normalized using the total number of genes found per MAG.

Cellulose-related CAZymes found within all MAGs in rice straw

Figure S4. Clustered MAGs recovered from rice straw samples and their normalized counts for CAZy families associated with lignocellulose-degrading enzymes.The organization of the tree relies on the phylogenetic information retrieved from GTDB-Tk2. The displayed MAGs are representatives from the generated clusters using ANI, and they were selected based on the following score: completeness – 0.5 x contamination.

Enzymatic activity of the different samples

Figure S5. Experimental protease, acid and alkaline phosphatase activities of rice straw (a) and soil (b) samples. The enzymatic activity was measured according to the guidelines established by Alef & Nannipieri (1995).

Physicochemical parameters of soil samples

Figure S6. Experimental measures of physicochemical parameters in soil samples including total Carbon, total Nitrogen, pH and phosphorus. The data were normalized using z-score to adjust the raw values to the same scale; the displayed numbers represent the mean of the quantified parameter per sample.